 function [imp_vect_relat, imp_vect_absol, ...
             F_stat_test1, p_values_test1, signif_test1, Rsquared_test1, Betas_test1, ...
             F_stat_test2, p_values_test2, signif_test2, Rsquared_test2, ...
             correl_instrument_shock, a_vect, estim_shocks] = zFC_IdentifExtInstruments5(Resid, Instruments, position_indicator, Sigmahat)

% -----------------------------------------------------%
% DESCRIPTION OF THE FUNCTION:
% this function identifies the structural VAR model using "Instruments" as
% external instruments 
%  
% INPUTS:
% 1) Resid: reduced form innovations of the underlying VAR
% 2) Instruments: set of one or more instruments to identify one underlying
% structural shock
% 3) position_indicator: sets the variable that will be used to define the
% normalization of the relative impulse vector and the absolute impulse
% vector
% 4) Sigmahat: variance-covariance matrix of the reduced form innovations
% of the VAR
% 
% OUTPUTS:
% 1) imp_vect_relat: relative impulse vector, normalized to 1 at equation
% "position_indicator"
% 2) imp_vect_absol: absolute impulse vector, where the scaling parameter
% is taken with positive value and is computed from the covariance
% restrictions of the VAR
% 3 to 7) F_stat_test1, p_values_test1, signif_test1, Rsquared_test1, Betas_test1: 
% results of the test on the strenght of the instrument(s), when regressing 
% residuals on instruments
% 8 to 11) F_stat_test2, p_values_test2, signif_test2, Rsquared_test2: same as above, but
% when regressing instrument(s) on residuals
% 12) correl_instrument_shock: correlation(s) of the instrument(s) with the
% underlying unobserved shock (respects normalization of the impulse vector)
% 13) a_vect: vector to comput structural shocks from reduced form shocks
% (respects normalization of the impulse vector)
% 14) estim_shocks: estimated structural shocks
% -----------------------------------------------------%
                                                                 
                                                                 

    % ============================ %
    % PRELIMINARY PART
    % ============================ %
    
    % Make sure variables enter as column vectors
    if size(Resid,1) < size(Resid,2) 
        Resid = Resid'; 
    end
    if size(Instruments,1) < size(Instruments,2) 
        Instruments = Instruments'; 
    end

    % Make sure they have the same length
    assert(size(Resid,1) == size(Instruments,1))

    % Create operational variables
    obs           = size(Resid,1);
    n_variab      = size(Resid,2);
    n_instruments = size(Instruments,2);


    % ============================ %
    % TESTS ON THE STRENGHT OF THE INSTRUMENT(S)
    % ============================ %

    % Using regressions of single residuals on instrument(s)
    F_stat_test1   = NaN*ones(n_variab,1);
    p_values_test1 = NaN*ones(n_variab,4);
    signif_test1   = NaN*ones(n_variab,4);
    Rsquared_test1 = NaN*ones(n_variab,1);
    Betas_test1    = NaN*ones(n_variab,1);

    for i = 1:n_variab;
        yyy      = Resid(:,i);
        xxx_rest = ones(obs,1);
        xxx = [ones(obs,1), Instruments];

        [F_test, Autocorr_test, White_test]  = zFC_Ftest_Autocorr_Homosk(yyy,xxx,xxx_rest);
        F_stat_test1(i)     = F_test(1);
        p_values_test1(i,1) = F_test(2);
        p_values_test1(i,2) = Autocorr_test(1,2);
        p_values_test1(i,3) = Autocorr_test(2,2);
        p_values_test1(i,4) = White_test(2);


        signif_test1(i,1) = zFC_significance_stars(F_test(2));
        signif_test1(i,2) = zFC_significance_stars(Autocorr_test(1,2));
        signif_test1(i,3) = zFC_significance_stars(Autocorr_test(2,2));
        signif_test1(i,4) = zFC_significance_stars(White_test(2));
        
        
        coeff = (xxx'*xxx)^(-1)*xxx'*yyy;        
        fitted = xxx*coeff;
        resid = yyy - fitted;
        yyy_demean = yyy - mean(yyy);
        Rsquared_test1(i) = 1 - resid'*resid/(yyy_demean'*yyy_demean);
        
        if n_instruments == 1
            Betas_test1(i) = coeff(2);
        end
    end


    % Using regressions of instrument(s) on all residuals
    F_stat_test2     = NaN*ones(n_instruments,1);
    p_values_test2   = NaN*ones(n_instruments,4);
    signif_test2     = NaN*ones(n_instruments,4);
    Rsquared_test2   = NaN*ones(n_instruments,1);
    
    for i = 1:n_instruments
        yyy      = Instruments(:,i);
        xxx_rest = ones(obs,1);
        xxx = [ones(obs,1), Resid];

        [F_test, Autocorr_test, White_test]  = zFC_Ftest_Autocorr_Homosk(yyy,xxx,xxx_rest);

        F_stat_test2(i)     = F_test(1);
        p_values_test2(i,1) = F_test(2);
        p_values_test2(i,2) = Autocorr_test(1,2);
        p_values_test2(i,3) = Autocorr_test(2,2);
        p_values_test2(i,4) = White_test(2);
        
        
        signif_test2(i,1) = zFC_significance_stars(F_test(2));
        signif_test2(i,2) = zFC_significance_stars(Autocorr_test(1,2));
        signif_test2(i,3) = zFC_significance_stars(Autocorr_test(2,2));
        signif_test2(i,4) = zFC_significance_stars(White_test(2));
    
        
        coeff = (xxx'*xxx)^(-1)*xxx'*yyy;        
        fitted = xxx*coeff;
        resid = yyy - fitted;
        yyy_demean = yyy - mean(yyy);
        Rsquared_test2(i) = 1 - resid'*resid/(yyy_demean'*yyy_demean);
    end
    
    
    % ============================ %
    % IDENTIFICATION OF THE MODEL
    % ============================ %
    
    % Check that the indicator variable corresponds to the residuals entering first, otherwise rearrange
    if position_indicator > 1
        % rearrange the order of the residuals
        step1 = [1:1:n_variab];
        step2 = step1 ~= position_indicator;
        step3 = step1(step2);
        Resid = [Resid(:,position_indicator), Resid(:,step3)];

        % rearrange the order of the Sigma matrix
        Sigmahat_step11 = Sigmahat(position_indicator,position_indicator);
        Sigmahat_step12 = Sigmahat(position_indicator,1:position_indicator-1);
        Sigmahat_step13 = Sigmahat(position_indicator,position_indicator+1:end);
        Sigmahat_step22 = Sigmahat(1:position_indicator-1, 1:position_indicator-1);
        Sigmahat_step23 = Sigmahat(1:position_indicator-1,position_indicator+1:end);
        Sigmahat_step33 = Sigmahat(position_indicator+1:end,position_indicator+1:end);
        
        Sigmahat = [Sigmahat_step11,  Sigmahat_step12,  Sigmahat_step13; ...
                    Sigmahat_step12', Sigmahat_step22,  Sigmahat_step23; ...
                    Sigmahat_step13', Sigmahat_step23', Sigmahat_step33];
        
        % create a vector used later to organize results according to the original order
        step4 = step3<position_indicator;
        step5 = step3(step4);
        step6 = step3>position_indicator;
        step7 = step3(step6);
        
        reorder = [step5'+1; 1; step7'];
    end
    
    % Estimate the relative impulse vector
    if n_instruments == 1 % Single regressions
        yyy      = Resid;
    %     xxx      = [Instruments, ones(obs,1)];
        xxx      = Instruments;
        beta_hat = (xxx'*xxx)^(-1)*xxx'*yyy;
        imp_vect_inconsist   = beta_hat(1,:)'; 

        % Compute relative impulse vector, normalizing the element that will enter the equation numbered "position_indicator"
        imp_vect_relat = imp_vect_inconsist/imp_vect_inconsist(1); % the policy variable has been reordered as first

    elseif n_instruments > 1 % 2SLS
        imp_vect_relat    = NaN*ones(n_variab,1);
        imp_vect_relat(1) = 1;
    
        yyy      = Resid(:,1);
%         xxx      = [ones(obs,1), Instruments];
        xxx      = Instruments;
        beta_hat = (xxx'*xxx)^(-1)*xxx'*yyy;
        fitted   = xxx*beta_hat;

        for j = 2:n_variab
            yyy = Resid(:,j);
%             xxx = [fitted, ones(obs,1)];
            xxx = fitted;
            gamma_hat = (xxx'*xxx)^(-1)*xxx'*yyy;
            imp_vect_relat(j) = gamma_hat(1);
        end
    end
    
    % Compute the scaling factor from the covariance restrictions and then the absolute impulse vector
    assert( size(Sigmahat,1) == length(imp_vect_relat)  ) 

    mu = imp_vect_relat(2:end);

    Sigma11 = Sigmahat(1,1);
    Sigma21 = Sigmahat(2:end,1);
    Sigma22 = Sigmahat(2:end,2:end);

    gamma = Sigma22 + Sigma11*mu*mu' - Sigma21*mu' - mu*Sigma21';
    b12_squared = (Sigma21 - Sigma11*mu)'*gamma^(-1)*(Sigma21 - Sigma11*mu);
    b11 = sqrt(Sigma11 - b12_squared); % take the positive solution for b11. If need to flip the sign of the entire b vector, change also the sign of the a vector and of the estimated correlation between instrument and structural shock of interest

    imp_vect_absol = imp_vect_relat*b11; 
 
    
    % Estimate the structural shocks 
    b21 = imp_vect_absol(2:end);
    b12primeB22inv = (Sigma21 - b11*b21)'*(Sigma22 - b21*b21')^(-1);
    a_vect_prime = (b11 - b12primeB22inv*b21)^(-1)*[1, -b12primeB22inv];
    a_vect = a_vect_prime';
    estim_shocks = Resid*a_vect;
    

    % Estimate the correlation of the instrument(s) with the unobserved structural shock
    correl_instrument_shock = NaN*ones(n_instruments,1); % one for each instrument given
    for i = 1:n_instruments
        E_r_m = (Resid'*Instruments(:,i))/length(Instruments(:,i)); % E(r_t * m_t). If position_indicator = 1, will have flipped order. phi is unaffected       
        step = E_r_m./imp_vect_absol;
        phi = step(unidrnd(n_variab)); % which element is taken is irrelevant, unless the first step also includes a constant and that we did not include, or unless we use more than one instrument it in the 2SLS. Otherwise, approximately the same
        correl_instrument_shock(i) = phi/std(Instruments(:,i)); % standard deviation of structural shock normalized to 1. If interpretation of the shocks requires flipping sign of imp_vect_absol, then the sign correl_instrument_shock , a_vect and estim_shocks needs to be changed
%         a_vect_bis = Sigmahat^(-1)*E_r_m/phi; % % if n_instruments=1, a_vect = a_vect
    end
    
    
    % Reorder the impulse vectors so that the policy indicator enters at the original position
    if position_indicator > 1
        imp_vect_relat     = imp_vect_relat(reorder);
        imp_vect_absol     = imp_vect_absol(reorder);
        a_vect             = a_vect(reorder);
    end  
    
end